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Abstract — Point  or  plot  data  are  commonly  available  for  mapping  forest  landscapes.  Be- 
cause such  data  are  sampled,  mapping  a  complete  coverage  usually  requires  some  type  of 
interpolation  between  plots.  SAS/GRAPH®  software  includes  the  G3GRID  procedure  for 
interpolating  or  smoothing  this  type  of  data  to  map  with  G3D  or  GCONTOUR  procedures. 
However,  the  smoothing  process  in  G3GRID  is  not  easily  controlled,  nor  can  it  be  used  to 
display  missing  data  within  rectangular  grid  maps.  These  shortcomings  motivated  devel- 
opment of  SAS  code  that  prepares  point  data  for  display  in  mapping  units.  This  code  links 
well  with  the  rest  of  the  SAS  system  to  allow  for  powerful,  easily  controlled  data  analysis 
within  mapping  units.  Examples  are  given  for  mapping  forest  vegetation  with  the  GMAP 
procedure. 
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Forest  landscapes  are  usually  described  using  sample  points  (plots)  lo- 
cated according  to  some  statistical  design.  In  order  to  map  these  data,  it 
is  often  desirable  to  interpolate  among  the  sample  plots  for  more  complete 
map  coverage.  SAS/GRAPH®  software  includes  the  G3GRID  procedure  for 
creating  a  rectangular  grid  of  interpolated  or  smoothed  values  from  ir- 
regularly spaced  observations  by  using  a  variety  of  spline  techniques.  These 
smoothed  data  can  be  mapped  with  G3D  or  GCONTOUR  procedures  to 
show  patterns  across  a  landscape. 

The  smoothing  process  of  G3GRID  is  not  easily  controlled,  nor  can  it  be 
used  to  display  missing  data  within  rectangular  grid  maps.  There  is  no 
control  for  allowable  distance  between  nearby  data  points  in  the  interpo- 
lation; only  the  number  of  nearby  data  points  can  be  controlled.  In  other 
words,  every  point  in  the  rectangular  grid  receives  a  value  based  on  its 
nearest  data  points,  regardless  of  the  distance  between  points.  Not  having 
some  distance  control  is  particularly  troublesome  for  mapping  fragmented 
landscapes  sampled  at  different  intensities.  For  these  situations,  it  is 
more  accurate  for  the  rectangular  grid  to  depict  nonsampled  areas  with 
missing  values  instead  of  assigning  values  estimated  from  distant  data 
points. 

The  shortcomings  in  G3GRID  procedure  motivated  development  of  sev- 
eral SAS  DATA  steps  that  summarize  point  data  into  mapping  units  for 
display.  This  paper  discusses  the  construction  of  mapping  units  and  illus- 
trates their  use  in  making  forest  vegetation  maps. 


This  paper  was  presented  at  SAS  Users'  Group  International  20,  Orlando,  FL,  April  2-5,  1995. 

David  C.  Chojnacky  is  a  Research  Forester  and  Mark  E.  Rubey  is  a  Computer  Programmer, 
Intermountain  Research  Station,  located  at  the  Forestry  Sciences  Laboratory,  Ogden,  UT  84401. 
e-mail:  /s=d.chojnacky/oul=s22102a@mhs-fswa.attmail.com. 
Telephone:  (801)  625-5402,  FAX:  (801)  625-5723. 
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Mapping  Units 


A  mapping  unit  is  the  fundamental  building  block  used  to  construct 
maps.  It  is  characterized  by  a  group  of  data  points  that  can  be  geographi- 
cally referenced.  Within  the  SAS  system,  a  mapping  unit  is  similar  to  a 
county  or  State  designation,  except  its  size  and  shape  is  determined  by 
the  spatial  scale  of  available  data  and  by  analysis  need.  Sample  data 
within  mapping  units  are  analyzed  with  standard  statistical  procedures 
to  obtain  summary  statistics  that  can  be  mapped  to  show  landscape 
patterns. 

To  illustrate  construction  and  display  of  mapping  units,  forest  vegeta- 
tion data  from  the  Kaibab  National  Forest  in  Arizona  are  used.  These 
data  include  a  mosaic  of  8,183  forest  stands  that  are  geographically  refer- 
enced at  each  stand's  center  using  the  Universal  Transverse  Mercator 
(UTM)  map  reference  system  (fig.  1).  A  stand  is  an  area  containing  simi- 
lar forest  vegetation  that  can  range  in  size  from  1  ha  to  over  100  ha. 
Stands  are  described  by  field  measurements  of  the  structure  and  composi- 
tion of  their  tree  species. 

Because  field  measurements  are  costly  and  are  made  over  long  time  pe- 
riods, it  is  common  to  have  situations — as  found  in  the  Kaibab  data  base 
(fig.  1) — where  all  contiguous  stands  have  been  identified  with  geographic 
coordinates  (from  maps)  but  not  all  stands  have  been  field  measured. 
Therefore,  it  is  important  that  the  mapping  techniques  for  forest  stand 
data  have  considerable  flexibility  for  dealing  with  missing  data. 

The  purpose  of  the  example  maps  is  to  describe  forest  conditions  around 
nest  sites  of  a  hawk  species  found  in  the  Kaibab  National  Forest.  In  Arizona 
and  New  Mexico,  the  northern  goshawk — a  forest  inhabitant — is  listed  by 
the  Forest  Service  as  a  "sensitive  species"  that  requires  special  manage- 
ment attention.  Maps  are  useful  for  examining  foraging  areas  for  gos- 
hawk home  ranges  scattered  throughout  forested  landscapes. 

Two  different  types  of  mapping  units  will  be  illustrated.  Kaibab  stand 
data  are  grouped  into  a  rectangular  grid  to  show  an  entire  landscape,  and 
stand  data  are  grouped  into  overlapping  circles  that  correspond  to  the  for- 
aging area  of  a  goshawk  pair  around  their  nest  site.  Only  key  concepts  for 
the  grid  and  circular  mapping  techniques  are  discussed  here,  but  more 
detailed  SAS  code  is  available  from  the  authors  upon  request. 

Rectangular  Mapping  Units 

In  a  DATA  step,  all  available  stand  data  are  placed  into  rectangles  (in 
this  case  squares)  of  a  fixed-size,  equally  spaced  across  the  landscape.  The 
center  or  centroid  of  a  stand  determines  to  which  mapping  unit  a  stand 
belongs.  The  size  of  the  mapping  unit,  whether  it  includes  few  or  many 
stands,  determines  the  map  scale.  Because  data  analysis  is  done  within 
each  mapping  unit,  there  is  flexibility  for  coping  with  data  problems  and 
land  use  patterns  that  might  be  encountered  in  mapping  a  natural  land- 
scape. For  example,  it  is  easy  to  isolate  and  treat  a  few  mapping  units  dif- 
ferently without  affecting  the  rest  of  the  map. 

Step  One — First  needed  to  make  a  map  are  the  four  X  and  Y  UTM- 
coordinates  (X0,  Y0,  Xe,  and  Ye)  that  border  the  entire  landscape.  Two 
nested  DO  statements  are  used  to  group  stand  data  into  mapping  units 
(MAPUNIT): 
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North  Kaibab  Ranger  District 


Figure  1 — The  North  Kaibab  Ranger  District  of  the  Kaibab  National 
Forest  in  Arizona  spans  almost  70  km  from  north  to  south,  and  it 
borders  the  Grand  Canyon  National  Park  on  the  south.  There  are 
8,183  forest  stands.  Most  of  the  inventoried  stands  include  the  high 
elevation  ponderosa  pine  (dry  site)  and  mixed  conifers  (moist  site). 
These  are  surrounded  on  the  east,  west,  and  north  by  very  dry  site 
pinyon-juniper  stands  that  mostly  have  not  yet  been  inventoried  (no 
field  data). 
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X  =  EASTING/1000;  *  a  stand's  UTM  X-coordinate  (1,000  m); 
Y  =  NORTHING/1000;  *  a  stand's  UTM  Y-coordinate  (1,000  m); 


XO  =  370; 
Xe  =  406; 
Y0  =  4020; 
Ye  =  4080; 


*  beginning  of  X-coordinate  range  (1,000  m); 

*  end  of  X-coordinate  range  (1,000  m); 

*  beginning  of  Y-coordinate  range  (1,000  m); 

*  end  of  Y-coordinate  range  (1,000  m); 


g=1.5; 


*  1,500-m  grid  size; 


Do  i  =  X0  to  Xe  By  g;    *  placing  X,Y  stand  data  into  mapping  units; 
Doj  =  Y0  to  Ye  By  g; 

If  (X  GE  i  and  X  LT  (i  +  g))  and 

(Y  GE  j  and  Y  LT  (j  +  g))  Then  MAPUNIT  =  Compress(i  I  I '-'  I  I  j); 


Step  Two— In  another  DATA  step,  the  center  or  centroid  of  each  map- 
ping unit  is  generated  for  all  possible  mapping  units  in  the  landscape  of 
interest: 

Do  i  =  X0  to  Xe  By  g; 
Doj  =  Y0  to  Ye  By  g; 

Xc  =  i  +  (g/2);  *  X-coordinate  for  MAPUNIT  center; 
Yc  =  j  +  (g/2);        *  Y-coordinate  for  MAPUNIT  center; 

MAPUNIT  =  Compress(i  I  I '-'  I  I  j);  *  for  identifying  MAPUNIT; 

If  (Xc  LT  Xe  and  Yc  LT  Ye)  Then  Output; 
End; 
End; 

This  facilitates  actual  mapping,  and  it  provides  a  way  to  identify  mapping 
units  that  may  have  no  field  data.  Results  from  this  step  are  merged  with 
step  one  results  (fig.  2). 

Step  Three — Analysis  is  done  within  mapping  units  by  using  SAS  sta- 
tistical procedures.  For  example,  frequencies,  totals,  means,  maximums, 
variances,  percentiles,  and  many  other  statistics  can  be  computed  to  sum- 
marize each  mapping  unit.  Because  the  stands  used  for  this  example  var- 
ied in  size  from  several  ha  to  over  100  ha,  it  was  useful  to  compute  map- 
ping unit  statistics  weighted  by  stand  size.  In  general,  use  of  the  SAS 
procedures  to  characterize  mapping  units  was  found  to  be  a  very  powerful 
analysis  option  that  is  lacking  in  many  GIS  packages. 

Step  Four — Data  produced  from  step  three  is  in  a  three-dimensional  for- 
mat that  can  be  mapped  with  G3D  or  GCONTOUR  procedures,  but  these 
procedures  include  additional  data  interpolation  among  mapping  units  that 
may  not  be  desirable.  On  the  other  hand,  the  GMAP  procedure  can  be  used  to 
produce  a  map  that  describes  mapping  units  exactly  as  they  were  computed. 
To  use  GMAP,  the  four  corners  of  each  MAPUNIT  are  output. 

h  =  g/2;  *  half  the  grid  width; 


End; 


End; 


X  =  Xc  -  h; 
Y  =  Yc-h; 
Output; 


*  lower  left  coordinate  of  MAPUNIT; 


X  =  Xc  +  h; 
Y  =  Yc  -  h; 
Output; 


lower  right  coordinate  of  MAPUNIT; 
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Figure  2 — Centers  or  centroids  of  mapping  units  for  a  1 ,500  m  grid 
constructed  across  the  North  Kaibab  Ranger  District. 


X  =  Xc  +  h; 

Y  =  Yc  +  h; 
Output; 

X  =  Xc  -  h; 

Y  =  Yc  +  h; 
Output; 


*  upper  right  coordinate  of  MAPUNIT; 


*  upper  left  coordinate  of  MAPUNIT; 

This  produces  a  file  analogous  to  the  SAS  supplied  county  or  State  data 
files  that  are  easily  mapped  by  the  GMAP  procedure. 
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Circular  Mapping  Units 


Example  Maps 


As  an  alternative  to  mapping  the  entire  landscape,  it  might  be  useful  to 
restrict  vegetation  mapping  to  special  areas  such  as  foraging  areas 
around  goshawk  nests.  For  a  goshawk  pair,  their  foraging  area  roughly 
corresponds  to  a  2,750  m  circle  around  their  nest  (Reynolds  and  others 
1992). 

Step  One — Forage  area  mapping  units  (MAPUNIT)  are  created  by  us- 
ing ARRAY  statements: 

X  =  EASTING/1000;  *  a  stand's  UTM  X-coordinate  (1,000  m); 
Y  =  NORTHING/1000;      *  a  stand's  UTM  Y-coordinate  (1,000  m); 

Array  F  F1-F71;  *  mapping  units  are  71  goshawk  foraging  areas; 

Array  FY  FY1-FY71;  *  Y-coordinate  of  MAPUNIT  (1,000  m); 
Array  FX  FX1-FX71;        *  X-coordinate  of  MAPUNIT  (1,000  m); 

Array  A_sq  A_sql-A_sq71;  *  triangle  height  for  MAPUNIT  radius; 
Array  B_sq  B_sql-B_sq71;  *  triangle  length  for  MAPUNIT  radius; 
Array  C  C1-C71;  *  triangle  hypotenuse  for  MAPUNIT  radius; 

R  =  2.75;  *  radius  of  goshawk  MAPUNIT  (1,000  m); 

Do  i  =  1  to  71;  *  assigning  stands  to  forage  areas  (MAPUNIT); 

A_sq{i}  =  Abs(Y-FY{i})*Abs(Y-FY{i}); 
B_sq{i}  =  Abs(X-FX{i})*Abs(X-FX{i}); 
C{i}  =  Sqrt(A_sq{i}  +  B_sq{i}); 

If  c{i}  LE  R  Then  MAPUNIT  =  F{i}; 
If  c{i}  LE  R  Then  Output; 
End; 

In  this  example,  it  is  possible  to  include  stands  in  more  than  one  mapping 
unit  when  the  distance  between  nests  is  smaller  than  the  diameter  of  the 
foraging  areas.  Application  of  this  logic  to  other  problems  requiring  very 
large  array  dimensions  could  be  complicated  by  slow  processing  or  AR- 
RAY statement  limitations,  but  it  works  fine  for  71  array  elements. 

Step  Two — As  in  step  three  for  the  rectangular  mapping  units,  many 
different  SAS  procedures  can  be  used  to  produce  statistics  describing  each 
mapping  unit. 

Step  Three — If  the  circular  mapping  units  are  simply  displayed  as 
circles,  overlap  is  a  problem  for  the  goshawk  nest  data.  But  Thiessen  poly- 
gons offer  a  convenient  way  to  divide  up  the  landscape  where  the  circular 
mapping  unit  data  can  be  displayed  in  a  polygon  around  each  nest.  A 
Thiessen  polygon  includes  all  the  area  closest  to  each  nest  relative  to  all 
the  other  nests.  Mathematically,  Thiessen  polygons  are  defined  by  the 
perpendicular  bisectors  of  lines  drawn  between  all  nests.  Because  the  SAS 
system  does  not  include  a  Thiessen  polygon  algorithm,  the  polygon  bound- 
aries were  created  using  the  ARC/INFO  geographic  information  system 
and  then  transported  into  SAS. 


For  producing  vegetation  maps,  the  forest  structure  within  each  map- 
ping unit  was  described  using  a  stand  density  index  (SDI)  (Husch  and 
others  1982;  McCarter  and  Long  1986).  This  variable  combines  average 
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tree  size  and  tree  density  of  each  stand  into  a  single  statistic  that  quanti- 
fies forest  vegetation.  Forest  vegetation  was  divided  into  (1)  dry-site  pon- 
derosa  pine  which  included  some  pinyon-juniper  and  other  woodland  spe- 
cies and  (2)  moist-site  mixed  conifer  and  aspen  species.  Overlaid  onto  the 
vegetation  maps  were  northern  goshawk  nest  sites  identified  by  the  aver- 
age number  of  fledglings  (0  to  3)  produced  each  year  (from  3  successive 
years  of  data). 

A  map  at  a  1,500  m  scale  using  rectangular  mapping  units  and  the 
GMAP  procedure  shows  the  highest  stand  density  indices  among  the 
moist-site  tree  species  (fig.  3).  It  is  interesting  that  the  most  productive 
goshawk  nests  (with  three  fledglings)  are  found  in  dry-site  forest  vegeta- 
tion with  moderate  or  low  vegetation  density  (med  or  low  SDI). 

Because  the  data  produced  for  the  rectangular  mapping  units  are  stan- 
dard three-dimensional  data  (fig.  2),  they  can  be  displayed  with  other 
three-dimensional  SAS  procedures.  For  example,  another  view  of  the 
figure  3  map  using  the  GCONTOUR  procedure  (fig.  4)  shows  few  nests  in 
the  highest  vegetation  density  (high  SDI),  and  it  shows  most  of  the  three- 
fledgling  nests  near  the  lower  boundary  of  the  moderate  density  vegeta- 
tion (med  SDI). 

It  is  easy  to  change  the  size  of  the  mapping  units  to  view  the  landscape 
at  different  scales.  For  example,  the  data  in  figure  3  mapped  at  a  5,000  m 
scale  omit  some  of  the  detail  but  show  similar  trends  (fig.  5). 

Because  a  goshawk  pair  will  center  its  hunting  activity  around  its  nest 
during  the  March  to  September  breeding  season,  it  is  useful  to  examine 
the  potential  foraging  area  around  each  nest  by  using  circular  mapping 
units  (fig.  6).  Differences  in  forest  structure  relate  to  differences  in  prey 
abundance.  It  is  interesting  that  none  of  the  highest  vegetation  densities 
(high  SDI)  on  moist-site  mapping  units  show  three  fledglings  produced. 

The  data  for  the  Kaibab  National  Forest  are  conveniently  within  one 
UTM  zone,  making  all  calculations  relatively  straightforward.  Had  the 
Forest  straddled  a  UTM  zone,  some  modifications  would  have  been  neces- 
sary. For  example,  in  another  study  the  same  mapping  technique  presented 
here  was  used  to  make  a  32  km  grid  map  for  all  of  Arizona  and  New  Mexico 
(fig.  7).  The  primary  modification  was  using  projected  radians  (obtained  from 
the  GPROJECT  procedure)  in  place  of  UTM  coordinates. 
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Figure  3— Average  stand  density  index  (SDI)  for  mapping  units 
constructed  by  using  1 ,500  m  grid  to  group  stands.  The  annual 
number  of  goshawk  fledglings  produced  at  each  nest  is  overlaid 
on  the  map. 
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Figure  4 — Average  stand  density  index  (SDI)  contour  lines 
corresponding  to  the  lower  ranges  for  the  low,  medium,  and  high 
SDI  classes.  The  annual  number  of  goshawk  fledglings  produced 
at  each  nest  is  overlaid  on  the  map. 
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low    SDI,     dry  site 

med  SDI  ,  dry  site  BBBi  med  SDI  ,  moist  site 
riicjri    SDI,     dry    site        HI   high    SDI,    moist  site 


Figure  5 — Average  stand  density  index  (SDI)  for  mapping  units 
constructed  by  using  a  5,000  m  grid  to  group  stands.  The  annual 
number  of  goshawk  fledglings  produced  at  each  nest  is  overlaid 
on  the  grid  map. 
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low    SDI,     dry    site        J5B8      med    SDI,    moist    site  k 
med    SDI,     dry    site        H   high    SDI,    moist  site 
high    SDI,     dry  site  

Figure  6— Average  stand  density  index  (SDI)  for  mapping  units 
constructed  by  grouping  stands  within  a  2,750  m  radius  around 
each  nest.  Because  the  circular  mapping  units  overlap,  Thiessen 
polygons  are  used  to  display  the  results.  The  annual  number  of 
goshawk  fledglings  produced  at  each  nest  is  overlaid  on  the  map. 
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Figure  7 — Annual  tree  growth  for  Arizona  and  New  Mexico's  pinyon- 
juniper  dryland  forest.  The  mapping  units  represent  average  growth 
computed  from  about  2  to  10  plots  per  32  km  grid  cell.  Because  X,Y- 
data  coordinates  spanned  three  UTM  zones,  the  map  was  made  by 
using  an  Albers'  equal-area  projection  (from  the  GPROJECT 
procedure)  of  the  UTM  coordinates  converted  (Chojnacky  and 
Tymcio  1987)  to  latitude  and  longitude  in  radians. 


Conclusions  

Because  the  G3GRID  procedure  has  limited  data  smoothing  capability, 
SAS  code  was  written  to  allow  flexible  mapping  of  point  data  with  GMAP 
and  other  SAS  procedures.  By  grouping  the  point  data  into  mapping  units 
as  a  first  step,  other  SAS  statistical  procedures  can  be  used  to  compute  a 
wide  range  of  statistics  for  map  display.  It  is  easy  to  map  at  different  spa- 
tial scales  and  display  data  in  different  shaped  mapping  units.  Example 
maps  demonstrated  displays  of  forest  structure  around  goshawk  nests. 
More  detailed  SAS  code  of  the  mapping  technique  is  available  from  the 
authors  upon  request. 
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